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Abstract. We provide a necessary and sufficient condition for existence of 
Gaussian cubature formulas. It consists of checking whether some overdeter- 
mined linear system has a solution and so complements Mysovskikh's theorem 
which requires computing common zeros of orthonormal polynomials. More- 
over, the size of the linear system shows that existence of a cubature formula 
imposes severe restrictions on the associated linear functional. For fixed preci- 
sion (or degree), the larger the number of variables the worse it gets. And for 
fixed number of variables, the larger the precision the worse it gets. Finally, 
we also provide an interpretation of the necessary and sufficient condition in 
terms of existence of a polynomial with very specific properties. 



1. Introduction 

In addition of being of self-interest, Gaussian quadrature formulas for one- 
dimensional integrals are particularly important because among all other quadra- 
ture formulas with same number of nodes, they have maximum precision. Indeed, 
when supported on p nodes, they are exact for all polynomials of degree at most 
2p — 1. Moreover, the nodes are exactly the zeros of the orthogonal polynomial of 
degree p and existence is guaranteed whenever the moment matrix of the associated 
linear functional is positive definite. 

When its multi-dimensional analogue exists (then called Gaussian cubature for- 
mula), it shares those two important properties of support and precision. That is, 
the nodes of a cubature formula of degree (or precision) 2p — 1 are the common 
zeros of all orthonormal polynomials of degree p. However, and in contrast to the 
one-dimensional case, its existence is not guaranteed even if the moment matrix of 
its associated linear functional is positive definite. 

To the best of our knowledge, the only necessary and sufficient condition avail- 
able is [4| Theorem 3.7.4] due to Mysovskikh, which states that an n-dimensional 
Gaussian cubature formula of degree 2m — 1 exists if and only if the orhonormal 
polynomials of degree m have exactly s m _i := common zeros. So this 

condition requires computing those zeros, or, equivalently, all joint eigenvalues of n 
associated Jacobi matrices; see e.g. [4j Theorem 3.6.2]. Equivalently, this amounts 
to check whether some n multiplication matrices of size [ n+m ^ ) commute pairwise; 
see [U Theorem 3.6.5]. 

On the other hand, in Cools et al. [1] one may find several lower bounds on the 
number of nodes required for a given specific precision; see e.g. [1] Theorem 9, 10, 
11,12, 13]. For more details on orthogonal polynomials and cubature formulas, the 
interested reader is referred to e.g. the excellent book of Dunkl and Xu [4] and the 
survey by Cools et al. pQ. 
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Contribution. This paper is concerned with existence of Gaussian cubature 
formulas associated with a Borel measure /i on R™, with all moments finite. Our 
contribution is to provide a necessary and sufficient condition for existence, different 
in spirit from (and complementary to) Mysovskikh's theorem which requires solving 
polynomial equations, or, equivalently, computing all joint eigenvalues of n Jacobi 
matrices. Namely, we prove that a Gaussian cubature formula or precision 2m — 1 
exists if and only if a certain overdetermined linear system has a solution. The 
coefficient matrix of this linear system comes from expressing the product P a Pp of 
any two orthonormal polynomials P a , Pp of degree to, in the basis of orthonormal 
polynomials (Pg), of degree at most 2m. In fact, in this expansion, only the constant 
coefficient and the coefficients of the degree-2m orthonormal polynomials matter. 
This linear system is always overdetermined and the larger n for fixed to (resp. the 
larger m for fixed n) the worse it gets. This shows that existence of a Gaussian 
cubature formula imposes severe restrictions on the associated linear functional. 

Finally, we also provide an interpretation of the necessary and sufficient condi- 
tion, namely that there exists a polynomial Q, which is a linear combination of 
orthonormal polynomials of degree 2to, and such that 

J P Q (x)P^(x)d M (x) = J P Q (x)P^(x)Q(x)d M (x), 

for all pairs (P a , Pg) of orthonormal polynomials of degree to.. 



2. Notation, definitions and preliminaries 

Let x = (xi, . . . , x n ) and let R[x] be the ring of real polynomials in the variables 
x and R[x]d its subset of polynomials of degree at most d, which is a vector space 
of dimension s<j := ( n ^ d ) ■ The number of monomials of degree exactly d is denoted 
by Td ■— A polynomial / G R[x]d can be identified with its vector of 

coefficients (f a ) = : f G R Sd , in the canonical basis (x Q ), a G N™. For any real 
symmetric matrix A, the notation A >; (resp. A >- 0) stands for A is positive 
semidefinite (resp. positive definite). 

Let us define the graded lexicographic order (Glex) a < g i j3 on N", which first 
creates a partial order by \a\ < \/3\, and then refines this to a total order by breaking 
ties when \a\ = \f3\ as would a dictionary with x\ — a, X2 — b, etc. For instance 
with n = 2, the ordering reads \,x\,x<2.,x\,x\x<2.,x\, 

For every integer to G N, let x m G R[x]*" m be the column vector of all monomials 
of degree to, with the < g i ordering. For instance, with n = 2 and k = 3, x 3 G M[x] 4 
and x 3 reads 

X = ix-^ , X^X2 , X\X^ , X2 ) 

With any sequence y = (y a ), a G N™, one may associate a linear function 
C y : R[x] ^ M by 

(2.1) / (= f«x a ) ^ £y(/) = E /« V°> f e R W- 

a a 

By a slight abuse of notation, and for any v G R[x] p denote by £ y (v) the vector 
(C y (vj)), j < p. Similarly, for any matrix V G R[x] pxr denote by -C y (V) the 
matrix (jC y (Vij)), i < p, j < r. If fi is a finite and positive Borel measure with 
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finite moments y = (j/q), ol G N n , then 

(2.2) C y (f) = f»V« = E f° I x "^(x) = / /d/i, 

t>£N" a£N" ^ ^ 

and is called a representing measure for y. 

Moment matrix. The moment matrix M<j(y) associated with a sequence y = 
(y a ), is the real symmetric matrix with rows and columns indexed by NJJ, and 
whose entry (a, (3) is just £ y (x Q+ ^) (= y a +p), for every a, [3 € NJJ- Alternatively, 
M^(y) = £y((l, x, . . . , x d ) T (l, x, . . . , x d )). Of course, if y has a representing mea- 
sure fj,, then M,i(y) y for all d G N, because 

(f,M d (y)f) = £ y (/ 2 ) = f fdfx>0, Vf er. 

However, the converse is not true in general, i.e., M<j(y) >; for all d, does not 
imply that y has representing measure. 

2.1. Orthonormal polynomials. Most of the material of this section is taken 
from Helton et al. [5]. Given [i and y = (y a ), a G N^, as in ()2.2|) . assume that 
Md(y) >- for every d. Then one may define the scalar product (•, -) y on R[x]d: 



(f,g) y := (f,M d (y)g> 



fgdfi, V/, S 6l[4 



(2.3) 



With e? G N fixed, one may also associate a unique family of polynomials {P a ) C 
M[x]d, a G N^, orthonormal with respect to fi, as follows: 

F Q G lin.spanjx' 3 : (3 < g i a} 

(P a ,Pp)y = S a =p, a,p GNS 
(F Q ,x' 3 ) y =0 if /3 < gl a 
(F Q ,x«) y > 0, oGNS, 

where J Q=l g denotes the usual Kronecker symbol. Existence and uniqueness of such 
a family is guaranteed by the Gram-Schmidt orhonormalization process following 
the 'Glex" order of monomials, and by positivity of the moment matrix Mj(y). 
See e.g. [1 Theorem 3.1.11, p. 68]. 

Computation. Computing the family of orthonormal polynomials is relatively 
easy once the moment matrix Md(y) is available. Suppose that one wants to 
compute p a for some a G NJJ . 

Build up the submatrix M^(y) obtained from Md(y) by keeping only those 
columns indexed by j3 < g i a and with rows indexed by a < g i a. Hence M^(y) 
has one row less than columns. Complete M^(y) with an additional last row 
described by (y) [a, 0\ — x* 3 , < g i a. Then, up to a normalizing constant, 
P a = dct M^(y). For instance with n = 2, d — 2 and a = (11), one has 



(2.4) 



x i y P(n)(x) = M dct 



2/oo 
2/io 
2/oi 

2/20 

1 



2/io 

2/20 

2/ii 

2/30 
X\ 



2/oi 
2/ii 

2/02 
2/21 
^2 



2/20 

2/30 

2/21 

2/40 
r 2 
•'T 



2/ii 

2/21 
2/12 
2/31 
X\X 2 



where the constant M is chosen so that J P, 2 u .d/i = 1; eee e.g. [HE]. 



4 



JEAN B. LASSERRE 



3. Gaussian cubature formula 

Given a finite Borel measure fj, on K C R™ with all moments finite, points 
x(i) £ R™, and weights 7* £ R, i = 1, . . . s, the linear functional / : — > R, 

(3.1) / ^ /(/):= £ 7 i/W*)), /G-F, 

is called a cubature formula of degree (or precision) p, if 1(f) = J K / d/i for all 
/ £ R[x] p ; the points (x(i)) are called the nodes. In other words, the approximation 
of the integral J fdfi by (|3.ip is exact for all polynomials of degree at most p. 

Cubature formulas are the multivariate analogues of quadrature formulas in the 
one dimensional case. By TchakalofFs theorem, given \x on K with finite moments, 
and d £ N, there exists a measure v finitely supported on at most points of K, 
and whose moments up to order at least d coincide with those of fi; see e.g. Putinar 
[5]. And so there exists a cubature formula of degree d. In the one-dimensional 
case, the celebrated Gaussian quadrature formula based on d points (called nodes) 
of K has precision 2d — 1 and positive weights. It exists as soon as the moment 
matrix of order d is positive definite, and all the nodes are zeros of the orthonormal 
polynomial of degree d. It is an important quadrature because among all quadrature 
formulas with same number of nodes, it is the one with maximum precision, 

When its multivariate analogue exists (then called a Gaussian cubature), the 
nodes are now the common zeros of all orthonormal polynomials of degree d, the 
precision is also 2c? — 1, and the weights are also positive. But its existence is 
not guaranteed even if all moment matrices are positive definite. For a detailed 
account on cubatures and orhogonal polynomials, the interested reader is referred 
to the very nice book of Dunkl and Xu [4], as well as the survey or Cools et al. Q]. 

3.1. Existence of Gaussian cubature. Mysovskikh's theorem 14, Theorem 3.7.4] 
states that a Gaussian cubature of degree m exists if and only if the orthonormal 
polynomials P m have exactly (" + '™ _1 ) common zeros. And so checking existence 
reduces to computing all joint eigenvalues of n so-called Jacobi matrices of size 
s m _i described in e.g. |4j p. 114]. We here provide an alternative criterion for 
existence which reduces to check whether a certain overdetermined linear system 
has a solution. 

Let /j be a finite Borel measure with support supp fj, = K C R™ , and with all 
moments y = (y a ), ol £ N™, finite. We will adopt the same notation in [4]. Let (P a ), 
a £ N n , be a system of orthonormal polynomials with respect to [a, and with each 
m £ N, let P m = [P H=m ] £ M[x] r -™ be the (column) vector (P Q<mi) , . . . , F a(mi . m) ) T , 
where (a' mi ') C NJ^ are all monomials of degree m, with cr mi ) < g i a^" 12 ^ ■ ■ ■ < g i 
a( mr ™\ For instance with n = 2, 

P2 = (P20, P11, Po2) T ■ 

For any sequence z = (z a ) denote by M<i(z) the moment matrix, with rows and 
columns index in N^, and with entries 

M d (z)[a,/3] = £ z (P a Pp), a,/3£N^, 

where C z : R[x] — >• R is the linear functional defined in (|2.ip . In other words, M(z) 
is the moment matrix expressed in the basis of orthonormal polynomials and can be 
deduced from Mj(z) by a linear transformation. Notice that if z = y then Mrf(y) 
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is the identity matrix I. 

Observe that for each 7,/3 £ NJ^, with \j\ = \j3\ = m, 



2m-l 



(3.2) P, Pp = A^ 2m) P + ]T A^ 2m) P i + A^' 2m) P 2m , 

for some row vectors A^j 2 ™" 1 G R rj , j = 0, 1, . . . , 2to. Let i m := r m (r m + l)/2. 

Theorem 3.1. Le£ A^ - 2 ™ 1 ' be the vector (A^ 2m) ) € K tm , and Zei A( 2m ' 2m ) G 
jUt m xr 2m ^ e ^ e ma i r j jX w hose rows are the vectors A^™' 2m \ defined in h3.2}) for 
7,/3 G N™ loitfc | 7 | = |/3| = m. 

There exists a Gaussian cubature formula of degree 2m — 1 if and only if the 
linear system 

(3.3) A(°' m) + A( 2m ' m ' u = 0, 
has a solution u G R r2m . 

Proof. Only if part. If there exists a Gaussian cubature (hence of degree 2m — 1), 
there is a finite Borel measure supported on the s m -i common zeros of P m ; see 
[U Theorem 3.7.4]. Let z = (z a ), a G N^, be the moments up to order 2m, of the 
Borel measure v. Since the cubature is exact for polynomials in M[x]2m-i ) 

z a = £ z (x«) = £ y (x Q ) = y a) Va G N£ m _ x . 

In addition, z/ being supported on the common zeros of P m , we also have 

^z(P 2 ) = J P a (x) 2 i/(dx) = 0, Va G N™ with |a| = m, 

that is, all diagonal elements of the positive semidefinite matrix £ z (P m P^) vanish, 
which in turn implies £ z (P m P^) = 0. Combining with (|3.2[) yields 

2m-l 

= £ Z (P 7 P,) = A 7 °/ m) + Y, A S 2m) ^(Pi) +A 7 7- 2m) £ z (P 2m ) 

3=1 =£ y (F;i)=0 

_ a (0,2m) « (2m, 2m) r , p s 



( A (0,2m) + A (2m,2m) ^(p^)) ^ 



which is (|3.3p with u = £ z (P2m)- 

If part. Conversely, assume that (|3.3[) holds for some vector u G W 2m . Consider 
the sequence z = (z a ), a G N^, with z a = y a for all a G NJm-i an d £z(P*) = w Q 
for all a G N^, with |a| = 2m. (Equivalently, £ z (P2m) = u.) And indeed such 
a sequence exists. It suffices to determine z a for all a G with |a| = 2m 

(equivalently, £ z (x 2m )). Recall the notation x fe = (x Q ), a G N n , with |a| = k, and 
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recall that for each d G N, (1, Pi, ... , Fd) is a basis of M[x](2, and so, 
(3.4) 



" 1 




1 








Pi 




Soi 


Si 





^2m 






Sl2m ' S( 2m -l)2m 


S2m 





" 1 




x 1 




^2m 



=s 



for some change of basis matrix S which is invertible and block-lower triangular. 
And so z solves 

u = s 2m £ z (x 2m ) + e/: z (i,x 1 ,... ,x 2m - 1 ) T , 



where 9 G R r2mXs(2m 1} is the matrix [S 02 m| • • • |S( 2m _i )2m ], and y 2m -i = (Va), 
\a\ < 2m - 1. And so £ z (x 2 ™)_= S^u - S^Qy^-i. 

Next, the moment matrix M m (z), which by dchnition reads 



M m _i(z) 



C z (¥ m ) £ z (P ro Pi) 3 
simplifies to 



■£ z (P m P^_i) | 



£ z (P m ) 

^z(PlP^) 



^z(P m PfJ 



(3.5) 



M m (z) 



1 P T ) 



because £ z (P;Pj) = £ y (P;Pj) = for all i ^ j with i + j < 2m - 1. But then we 



£ Z (P 7 P,) = A<°/ m > 



also have £ z (P TO P, n ) = because from 

2m— l 

]T A^ 2m) £ Z (P,) +A( 2 ; i ' 2m) £ z (P 2m ) 

= ( A ( o ' 2m )+A( 2m ^u)( 7 ,/3)=0, 

for all 7, p e N™ with | 7 | = \/3\ = m. 

And so M m (z) ^ 0, and rankM m (z) = rankM m _i(z). By the fiat extension 
theorem of Curto and Fialkow, z is the moment sequence of a finite Borel measure 
ijj on M™, supported on rankM m (z) = s m _i distinct points x(fc), k = 1, . . . , s m _i; 
see e.g. Curto and Fialkow [H|3], Lasserre Theorem 3.7], or Laurent [6]. That 
is, denoting by S x the Dirac measure at x € R", 

S771-1 

fc=l 

for some strictly positive weights jk ■ From = £ z (P m P^) we obtain 

= £ Z (P 2 ) = f P Q (x) 2 ^(x), Va e N™, with |a| = m, 
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which proves that each x(fc) is a common zero of P m , k = 1, . . . , s m -i. But by [4j 
Corollary 3.6.4], P m has at most s m _i common zeros. Therefore, P m has exactly 
s m _i common zeros and ip is supported on aZZ common zeros of P m . By [3] Theorem 
3.7.4], admits a Gaussian cubature formula of degree 2m— 1. Indeed, since z Q = y a 
for all \a\ < 2m — 1, then for every / G R[x]2 m -i, 

f /(x)d/*(x) = £ y (/) = £,(/) = f fdi>=Y j 7fc/(x(fc)), 

with 7^ > for every k. □ 



By [1] Theorem 7], we also know that the x(fc)'s belong to the interior of the 
convex hull conv K of K. 

Observe that existence of a solution u G M 1 " 2 " 1 for the linear system (|3.3[) imposes 
drastic conditions on the linear functional C y because p. 31) is always overdeter- 
mined. Moreover, for fixed m (resp. for fixed n) the larger n (resp. the larger m) 
the worse it gets. This is because t m (= r m (r m + l)/2) increases faster than r 2m . 



Interpretation. We next provide an explicit expression of the vector A^ 0,2 ™-' as 
well as the matrix _\( 2m ' 2m \ which permits to obtain an interpretation of the 
condition (|3.3[) . 

For 7, (3 G N™, applying the linear functional C y to both sides of (|3.2|) yields 

(3.6) <5 7=/3 = £ y (P 7 P^) = A(°/ m) , V| 7 | = |/3| = m. 

because £ y (Pfc) = for all k. This provides an explicit expression of the vector 
A( 0,2m '. Similarly, multiplying both sides of (|3.2j) by P K with \k\ = 2m, and 
applying again £ y , yields 

(3.7) C y (P,P p P K ) = £ A^ 2m )(a)£ y (F Q P K ) = A^' 2 "^), 

|a|— 2m 

because £ y (P Q P K ) = 8 a = K . This provides an explicit expression for all the entries 
of the matrix A( 2m,2m ). And we obtain: 

Corollary 3.2. There exists a Gaussian cubature of degree m if and only if there 
exists a polynomial Q — u T P2 TO , for some u G M. r2m , such that 

(3.8) CyfaPp) = £y(PyPpQ), 

for all /3,7 G N"j wzi/i I7I, |/3| = m. Equivalently, i3.8\) reads 

J P 7 (x)P (3 (x)^(x) = /" P 7 (x)P /3 (x)Q(x)d M (x), 
/or all /3,7 G NJ^ wrai/i I7I, |/3| = m 7 or m matrix form: 

1= /p,„(x)P m (x) T Q(x)^(x), 
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Proof. With u as in (JUSJ) and using (|3ll)-([3J]), for all /3, 7 G N™ with = m, 

one obtains 

C y (P,P p ) = A^ 2 " 1 ' = (A(»)u) %/3 

= £ y P 7 P^ ( ]T Uq p q ) 

V aGNJ m ,|a|=2m 

= Cy (P 7 P /9 (U T P 2m )) , 

which is (|3.8[) with Q = u T P 2m . Conversely, if (JSTSJ) holds with for some Q = u T P 2m 
then obviously u is a solution of (I3.3[) . □ 



Final remark. One may obtain more information about the polynomial Q of 
Corollary 13.21 From the proof of Theorem 13.11 for any solution u £ W 2m of (|3.3p , 
one has u = £ z (P 2m ) where z is the moment vector of the measure supported on 
the s m -i nodes of the Gaussian cubature. That is, 



u = £ z (P 2m ) - 7fcF 2m (x(fc)), 



fc=i 



for some positive weights 7&. On the other hand, the polynomial Q of Corollary 
is of the form u T P 2m where u solves (|3.3[) . Therefore, 



x 1— ► Q(x) = u T P 2m (x) = 7fcP 2 m(x(fc)) T P 2m (x). 

fc=l 

Hence Q satisfies 

£ y (P Q g) = u T £ y (P Q P 2m ) = 0, V\a\ < 2m 
(in particular J Qd/i = 0), and 

S m -1 

C y {P a Q) = u a = 7fc P Q (x(fc)), V|a| = 2m. 
fe=i 
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